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, Context. Synthesis models predict the integrated properties of stellar populations. Several problems exist in this field, mostly related to the fact 
■ that integrated properties are distributed. To date, this aspect has been either ignored (as in standard synthesis models, which are inherently 
f*"^ ' deterministic) or interpreted phenomenologically (as in Monte Carlo simulations, which describe distributed properties rather than explain 
! them). 

' Aims. This paper presents a method of population synthesis that accounts for the distributed nature of stellar properties. 
~"f-3 Methods. We approach population synthesis as a problem in probability theory, in which stellar luminosities are random variables extracted 
Q ■ from the stellar luminosity distribution function (sLDF). 

Results. With standard distribution theory, we derive the population LDF (pLDF) for clusters of any size from the sLDF, obtaining the scale 
relations that link the sLDF to the pLDF. We recover the predictions of standard synthesis models, which are shown to compute the mean 
of the luminosity function. We provide diagnostic diagrams and a simplified recipe for testing the statistical richness of observed clusters, 
thereby assessing whether standard synthesis models can be safely used or a statistical treatment is mandatory. We also recover the predictions 
of Monte Carlo simulations, with the additional bonus of being able to interpret them in mathematical and physical terms. We give examples 
' of problems that can be addressed through our probabilistic formalism: calibrating the SBF method, determining the luminosity function 
k>( , of globular clusters, comparing different isochrone sets, tracing the sLDF by means of resolved data, including fuzzy stellar properties in 
j_j . population synthesis, among others. Additionally, the algorithmic nature of our method makes it suitable for developing analysis tools for the 
' Virtual Observatory. 

Conclusions. Though still under development, ours is a powerful approach to population synthesis. In an era of resolved observations and 
pipelined analyses of large surveys, this paper is offered as a signpost in the field of stellar populations. 

Key words. Clusters - Galaxies: star clusters - Galaxies: stellar content - Hertzsprung-Russell (HR) and C-M diagrams - Methods: data 
analysis 
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1. Introduction and motivation 

The study of stellar populations is one of the most fecund topics 
in today's astronomy. Understanding the properties of stellar 
populations is a key element in the solution of a host of funda- 
mental problems, such as the calibration of distances to extra- 
galactic objects, the age determination of clusters and galaxies 
through color fitting, the characterization of the star-formation 
history of composite populations, the modeling of the chemical 
evolution of galaxies, and several more. When tackling these 
problems, the interaction between theory and observations goes 
both ways: one may want to predict the properties of a stellar 
population with given properties, or to recover the basic proper- 
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ties of an observed population from the available observables. 
In either case, to accomplish the task one needs feasible theo- 
retical models that serve as reliable diagnostic tools. 

A traditional approach to building such diagnostic tools has 
been the computation of synthesis models. Synthesis models 
allow one to predict the features and the evolution of stellar 
populations in a highly structured way, one that is apt for rou- 
tine analysis and quantitative assessments. For example, syn- 
thesis models can be used, in combination with other methods, 
for the determination of stellar population properties of large 
samples of galaxies in a reasonable time; e.g. the analysis of 
50362 galaxies of the Sloan Digit al Sky Survey based on thei r 
integrated properties performed bv lCid Fernandes et al.l J2005I) . 

However, in recent years there has been a growing aware- 
ness that synthesis modeling also suffers from severe limita- 
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tions. In some cases, we know a priori that standard mod- 
els cannot be applied to a given system, because the prop- 
erties of the system fall outside the hypothesis space of the 
models; this is the case, for example, of un dersampled popu- 
lations JChiosi et alJll988tlCervino & Valls-Gabaudll2003l) . In 
other cases, we observe a mismatch between the properties of 
the system inferred from observations of their individual com- 
ponents and those derived from their integrated properties an- 
alyzed with synthesis models: e.g. the discrepancy between 
the age determined from the color-magnitude diagram and the 
spectroscopic ag e in NGC 588 Jjamet et alJl2004 . or the IMF 
slope inferred hv lPelleritilll2004 in undersampled giant H n re- 
gions. In these cases, we are facing a crisis in the explanatory 
power of synthesis models. 

Previous attempts to solve this crisis and to understand the 
limitations of synthesis models have repeatedly pointed at the 
necessity of including statistical considerations in the analy- 
sis. According to this view, the predictive power breaks down 
because the traditional synthesis model is essentially a deter- 
ministic tool, whereas the nature of the problem is inherently 
stochastic. The clearest example of stochasticity is the mass 
spectrum of stellar populations, in which fluctuations (possibly 
random, although not necessarily so) in the number of stars of 
each type appear around the mean expected number. Until now, 
the efforts to take stochasticity into account in the modeling of 
stellar populations have moved in essentiall y two directions 
the use of Monte Carlo techniques (e.g . Santos & Frogel 1997 
Cj^ino J1 I :1 uridjajia^ iBru zual 2002; Girardi 

2002; Cantiell o et al.ll2003l among others) and the r einterpre- 
tation of traditional models in statistical terms (e.g. Buzzoni 
ll989HCervino et alJl2002bl iGonzalez et alJl2004 . Both meth- 
ods have proved able to solve some of the problems, or at 
least point toward possible solutions. However, they suffer from 
practical difficulties. Monte Carlo simulations are expensive (in 
terms of disk space, CPU time, and human time required to an- 
alyze the results), while, to date, the statistical reinterpretation 
of standard models has only served to establish limits to the use 
of synthesis model s for clusters with moderately under sampled 
populations JCervino & Luridianall2004[ ICervino et al.l Eo03). 
i.e. clusters with total initial masses on the order of 10 5 M . 

This limitation brings about a serious impasse for the study 
of stellar populations by means of their integrated light, since 
the properties of clusters with lower masses cannot be reliably 
obtained. This class includes, for example, all of the clusters in 
our Galaxy (including globular clusters), clusters in the Large 
Magellanic Cloud (LMC), as well as many clusters in external 
galaxies. F or example, many of the clusters in the 'Antennae' 
studied by Zhang & Fal] II 19991) have masses lower than this 
limit: in Sect. 18.31 we will show that explicit consideration of 
stochastic effects could alter the conclusions that are drawn 
on the cluster luminosity function based on these clusters. 
Furthermore, these limitations will become even more dramatic 
with the development of Virtual Observatory (VO) technolo- 
gies, which will make the automatic analysis of large amounts 
of data possible. It is therefore mandatory to adapt evolution- 
ary synthesis models to the present needs, so that they can be 
applied to observations. 



In the present paper we introduce a theoretical formal- 
ism for the probabilistic analysis of single stellar populations 
(SSPs). Our formalism yields results that are as accurate as 
those of large Monte Carlo simulations, but it bypasses the need 
to perform these simulations: that is, the method is both accu- 
rate and economic. By means of our formalism, synthesis mod- 
els can be applied to clusters of any size and the confidence 
intervals of the results can be evaluated easily. This makes it 
possible to estimate the relative weights of different bands for 
the realistic application of goodness-of-fit criteria like the x 2 
test. Finally, the algorithmic nature of our method makes it fea- 
sible for implementation in the VO environment. 

This paper is the fourth in a series 
llCervifio. Luridiana & Castande^l200(i ICervino et alJl2002rt 
ICervino et alJl2001al) dealing with the statistical analysis of 
stellar populations. Although it only deals with SSPs, a fifth 
paper, in preparation, will be devoted to the extension of this 
formalis m to any sta r form ation history scenario. Finally, in 
ICervino & Luridiana! (|2005) we give an extensive review of 
the uncertainties affecting the results of synthesis models. In 
addi tion to th e se p apers, we also recommend the work by 
Gilfanov et al. (2004): Statistical properties of the combined 
emission of a population of discrete sources: astrophysical 
implications. This paper, although not directly focused on 
synthesis models, suggests an alternative point of view of the 
problem. In some aspects, it has inspired the present paper. 

In this work we restate the problem of stellar population 
synthesis from a new perspective, the one of luminosity dis- 
tributions, which is a powerful and elegant way to understand 
stellar populations. We provide several examples of applica- 
tions to illustrate such power and indicate potential areas of 
future development. The starting point of the new formalism is 
the definition of the stellar (i.e. individual) luminosity distribu- 
tion function (LDF) and of its relation to the central problem of 
population synthesis (Sect.|5Jl. Sect. [3] describes the two main 
variants of synthesis models by means of which the problem 
has traditionally been tackled: Monte Carlo and standard mod- 
els. Next, we take the reader on a journey through the no man's 
land of population synthesis' pitfalls (Sect. |4j, and show how 
these are faced by existing techniques (Sect. [3J. This account 
is necessary to introduce, in Sect. [6] our suggested solution 
and show its comparative power. Finally, in Sect. [7] we give a 
few examples of applications of our formalism to current prob- 
lems. Future developments of the present work are described in 
Sect. 03 and our main conclusions are summarized in Sect. [5] 



2. Overview of the problem 

This section starts by introducing a few basic definitions. An 
exhaustive summary of the notation used and its rationale is 
offered in AppendixlAl 

The general problem approached by synthesis models is the 
computation of the luminosity 1 L tot emitted by an ensemble of 
Mot sources - a stellar population. From a theoretical point of 
view, this problem can be characterized in three basic ways. 



1 Throughout the paper, 'luminosity' will be a generic label for the 
stellar emission in any given band. 
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If the luminosities Z, of the individual sources are known, 
the total luminosity L tot is obtained trivially as the sum of all 
the /, values: 



' ' tot 



(1) 



1=1 



This circumstance is not very frequent. Its most common exam- 
ples are Monte Carlo synthetic clusters (Sect.01 in the theoreti- 
cal domain and resolved stellar populations in the observational 
one. 

In the more usual situation in which the luminosities of 
individual objects cannot be specified on an individual basis, 
a different approach must be adopted: in this case a luminos- 
ity distribution function (LDF) tfihiO is assumed that describes 
the distribution of luminosity values in a generic ensemble. 
Traditionally, the LDF has been seen as the asymptotic limit 
of a differential count of stars: 



(p h (£) = lim 



(— lim ^ 

A M ot Af ^° 



(2) 



where AN is the number of stars in a luminosity bin of width 
Af, and N tol the total number of stars. The integral of the LDF 
is normalized to 1 : 



f 

Jo 



<p h (€)d€=l. 



(3) 



The integrated luminosity of an ensemble is traditionally ob- 
tained by means of the expression: 



Jo 



(4) 



The bottom line of the present work is that this approach is 
conceptually wrong and operationally sterile. The crucial point 
where we part company with this approach is the definition and 
interpretation of the LDF: to us, fhiO is a probability density 
function (PDF) from which the luminosity of an actual sys- 
tem is drawn. If the system is an individual star, its PDF is the 
stellar LDF (sLDF). If it is a stellar population, its PDF is the 
population LDF (pLDF). Hereafter, we will indicate the sLDF 
with the symbol (Pl(() and the pLDF with the symbol tfiL ml (-0- 

To avoid confusion, note that the pLDF is not the same as 
the cluster or galaxy luminosity function (LF) as commonly de- 
fined, because galaxy LFs include the effect of a spread in ages 
and number of stars (or, equivalently, ages and mass), whereas 
the pLDF defined in this paper represents the luminosity distri- 
bution of ensembles of stars with the same physical parameters 
(age and total number of stars). 

In this paper, we are concerned with the properties of the 
sLDF and its relation to the pLDF. In particular, we will show 
that the total luminosity of a stellar population, is a dis- 
tributed quantity, whose mean value M' x is given by: 



Mi = {£) = N« 



Jo 



l<p L {€)d€. 



(5) 



The integral on the right-hand side of this equation is the mean 
value of the sLDF, p\, that is: 



M\ = A^tot/ii- 



(6) 



Although Eqs. |4] and [5] yield similar expressions, it is im- 
portant to distinguish between the two approaches. The his- 
torical origin of our sLDF lies in star counts, but it is a step 
further from them, in the same sense that the frequentist and 
the objectivist definitions of probability differ: the frequentist 
definition depends on the realization of trials, while the ob- 
jectivist definition assumes that the probability properties are 
built-in. The implications of either approach will be discussed 
at length throughout the paper. For the time being, note that 
Eq. ^ involves a discrete sum, while Eqs. |4] and [5] involve an 
integral, whose numerical solution requires binning the inde- 
pendent variable. In Sect. 14. 21 we discuss the consequences of 
binning distribution functions for numerical applications. 

The main goal of synthesis models is to obtain the luminos- 
ity of a model stellar population, either by direct count (Eq.^ 
or by an integral including the sLDF (Eqs.|^and|5|i. In the fol- 
lowing, we will see how this can be carried out in practice, un- 
der the assumption that there has been a star-forming episode 
in which all the stars have formed simultaneously; this is the 
scenario assumed by SSP models. 

Because stellar luminosities evolve with time, the sLDF is 
a function of the star's age. Since the luminosity evolution of 
a star depends on its initial mass, we can express the time de- 
pendence of the sLDF explicitly by writing the sLDF as a func- 
tion of two other functions: the isochrone Km; t) and the initial 
mass function (IMF) ify^ni). The isochrone i(m; t) gives the lu- 
minosity of a star as a function of its initial mass m at a given 
value of the age t. The IMF gives the probability distribution of 
initial stellar masses. The IMF has a status similar to that of the 
sLDF, in that it can be either interpreted as the result of a star 
count: 



/ 1 AN\ 
Wuim) — lim lim 1 



Am / ' 



(7) 



or, in probabilistic terms, as the probability density for a star 
of being born with mass m: for reasons similar to those dis- 
cussed above, we support the second interpretation (which, by 
the way, implies that the IMF should better be called Initial 
Mass Probability Density Function). According to its defini- 
tion, <fiM(m) fulfills the normalization condition: 



Jo 



(p M (m)dm = 1. 



(8) 



Summing up, the sLDF can be rewritten in terms of the IMF 
and the isochrone as follows: 



ip L ({; t) = <p M (m) x 



dl(m\ t) 
dm 



(9) 



Eq.|5]can now be used to rewrite the mean value of the sLDF 
in terms of the isochrone and the IMF 2 : 



p x {t) - I Km, t) <p M (m) | — : | — : dm 



dm 



dm 



((m; i) (fmifn) dm, 



(10) 



2 The isochrone is not monotonic, so that the integral limits of Eq.|5| 
do not correspond to those of Eg. 1101 
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where the integration variable in Eq. [5]has been changed from 
I to m, and m low , m up are the lower and upper mass limits re- 
spectively of the integration domain. 

Solving this integral is the main task of stellar population 
synthesis modeling. In the following section we will describe 
the main types of synthesis models and how they perform this 
integral. 

3. Basic strategies 

Once the physical problem is framed, we must translate it into 
actual computations, a task carried out by evolutionary syn- 
thesis codes. These come in two basic flavors, standard and 
Monte Carlo. Standard simulations are models in which the ini- 
tial mass mixture is analytically described by the IMF, whereas 
in Monte Carlo simulations the ensemble of stars is selected 
by random sampling the mass of each star to be included in 
the population, and the IMF is used as the weighting function. 
Therefore, the mass distribution obtained with a standard sim- 
ulation is univocally determined by the population's parame- 
ters, while in Monte Carlo simulations it is not; additionally, 
the standard distribution is smoother than the Monte Carlo one. 
Both methods, however, operate on a binned mass spectrum, 
due to the limited resolution of numerical computation and to 
the discreteness of the available stellar models. This fact has 
important consequences, which will be discussed in Sect. 14.21 
Using either of these two approaches, evolutionary synthe- 
sis codes aim to characterize the integrated emission properties 
of an ensemble of stars as a function of its physical parameters, 
such as the age and number of the individual stars of the en- 
semble 3 . Standard codes perform this task by carrying out the 
integration of Eq.^ll As described in Sect. [2] the result can be 
interpreted in two alternative ways: either a deterministic one - 
yUj is the sum of the luminosities of all the stars included in the 
ensemble modeled, normalized to one star -, or a probabilis- 
tic one - yu'j is the mean value of the sLDF Although the two 
interpretations may seem close at first sight, they are funda- 
mentally different, both from a conceptual and from a practical 
perspective: this point will be furthered in Sect. [6] We call these 
two alternative interpretations of standard synthesis models the 
deterministic and the statistical one. Note that these labels do 
not identify different classes of models, but rather different in- 
terpretations within the same class of models - standard mod- 
els. In practice, some codes do not explore this interpretation, 
while others acknowledge the distributed nature of luminosity 
and compute, in addition to the mean luminosity, the variance 
of the distribution. In either case, fx^ can eventually be scaled 
to the size of the ensemble by multiplying by N to t\ this property 
will be formally demonstrated in Sect. 16.21 

3 Strictly speaking, not all the luminosity sources contributing to 
the integrated luminosity of a stellar population need be stars, as 
they include, e.g., accretion disks or thermalization of kinetic energy. 
Accordingly, one should in principle use the more general expres- 
sion 'luminosity sources' rather than 'stars' to avoid loss of gener- 
ality. This distinction, however, is not relevant for our scope here, 
since in our treatment any luminosity source can be reconducted to 
a star. Therefore, in the following we will use 'stars' and 'luminosity 
sources' as synonyms. 



Note that many codes do not use tfMQri) in Eq.^3 but rather 
a proportional function tp' M (m) = const (pu{m) normalized in 
such a way that: 

m <p' M (m)dm — 1 . (11) 

o 

Since mtpM(m)dm is the mean mass value (m) of the IMF, 
using <p' M (m) instead of (pyi{m) in Eq.^5|yields the mean lumi- 
nosity of one star divided by (m): in this case, the mean total 
luminosity is found by multiplying by the total mass of the en- 
semble M tol instead of N tot . 

As for Monte Carlo models, their task is essentially the 
computation of the sum in Eq.^ Each time a simulation is run, 
a particular realization is drawn from the underlying distribu- 
tion. The result of the simulation is the integrated luminosity X. 
of that particular realization. If many Monte Carlo simulations 
are available for a fixed set of input parameters, an estimate for 
M' x can be obtained as the mean of the _£ values. This implies 
that the results of Monte Carlo simulations depend not only on 
the underlying luminosity distribution, but also on the number 
of simulations used to sample such distribution. If the set of 
simulations is sufficiently large, the Monte Carlo method has 
also the potential to provide the actual distribution function of 
the luminosities of the ensemble. Hence, an important draw- 
back of the method is that it is intrinsically expensive, because 
the accuracy of the results increases with the number of simu- 
lations. 

In summary, the goal of determining the luminosity of stel- 
lar populations reduces to the computation of a sum (Eq.^ or 
an integral (Eq.llOli. However straightforward this may seem, 
this computation is hindered in the practice by several intrinsic 
features of the problem: these will be the topic of next section. 

4. Pitfalls in the handling of the LDF 

Our previous discussion has taken place on an abstract level. 
In the practice of synthesis modeling it is necessary to translate 
the concepts discussed above into specific prescriptions for the 
handling of equations, and deal with the limitations imposed by 
the input ingredients, which have a finite resolution in the pa- 
rameter space. In this section we will discuss a specific aspect 
of this task, namely the difficulties inherent to the determina- 
tion of the mean value of the sLDF. As a first step, let us revisit 
a few well-known results in terms of the sLDF. 

4.1. Domain limits and average properties 

In synthesis modeling it is a well established fact that the in- 
tegrated luminosity of the ensemble is dominated by the most 
massive stars in the cluster. In the following we will illustrate 
this point by means of three simplified scenarios that are repre- 
sentative of real sLDFs. These are shown graphically in Fig.^ 
In each panel, the position of the mean fi[ and the region corre- 
sponding to yu'j + 1 cr are explicitly marked. 

Case 1 : Main sequence luminosity function In the first case 
let us assume that all the stars are in the main sequence (MS), 



M. Cervino, V. Luridiana, and N. Cervino-Luridiana: Probabilistic synthesis models 



5 



a = 2.35; (3 = 3 
,,,-8.5 ,, 



/ij ° = a = 2700 £ Q 
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10 J 0.01 



10 100 1000 10^ 10 = 10° 
Luminosity \£q\ 



Fig. 1. Schematic representation of the of the sLDF for three 
different cases: main sequence only (top), main sequence plus 
a constant post main sequence (middle), and main sequence 
plus a bumpy post main sequence (bottom). The position of 
the mean (//j) and the region corresponding to jj\ + 1 <j are 
explicitly marked. 

and that t oc mP. Assuming a power-law IMF with index —a, 
the sLDF can be expressed as: 



a 1 1 

w h oc£ ■ - ( = - ( . 

^ fi fi 



(12) 



(Sal peter|ll955h . If the IMF departs from a simple power law, 
the conclusions are less straigh tforward. For exampl e, if we 
consider the log-normal IMF by Mille r & Scald d 19791) and ap- 
proximate it as a power-l aw series, a val ue of a > 4 is reached 
near 100 M (see Fig. 1 in Kroup abOOll) . In this case, the dom- 
inant stars will not be the most massive ones. 



Case 2: Luminosity function with a constant post main- 
sequence contribution Let us go a bit further and add a post- 
main sequence (PMS) contribution to the sLDF. The sLDF can 
now be divided in two regimes, corresponding to the MS and to 
the PMS respectively. As a first approximation, we will assume 
that the PMS phase gives a constant contribution to the sLDF, 
implying that any point of the PMS portion of the isochrone 
is equally probable. The sLDF can be described by the expres- 



<Pl = 



l-a l-a 

p~T~ _ f ~T~ 



(15) 



A l max 

~ c„ r „-e- 



TO 



TO 



if ( € {(jo, Anax), 



where /'to is the turn-off luminosity. The absolute value of the 
PMS contribution is chosen so as to yield the same coefficient 
A for the two regimes, while preserving the overall normaliza- 
tion; that is, we are forcing the PMS to have a fixed weight 
with respect to the MS phase. This choice allows us to keep 
the complexity of the expressions to a minimum, but it has no 
influence on the general conclusions, which would be reached 
even if we dropped this assumption. 



If we further assume that ( T q 
of the sLDF is: 



» I 

mm 



, the mean value 



1 +p-a 



P 

TO 



2(1 -a) 



+ fro). (16) 



That is, the mean depends on both £jq and ( mix . Note that, 
since in the PMS phase the relation between the mass and the 
luminosity is not simple anymore, the above result cannot be 
easily expressed as a function of the initial mass. As a trivial 
example, if the age t is larger than the lifetime of a star of initial 
mass m up , then £(m ap , t) = 0. 



The constant factor can be determined by imposing the normal- 
ization condition: 



1 



1 



P p — r 0° 

'-max 



l-a-0 



A 1---0 



The mean value of the luminosity is then: 



Pi 



- 

— dt 



1+/3- 



i+0-' 
i 



■ C 



(13) 



(14) 



Case 3: sLDF with a bumpy post main-sequence contribu- 
tion As a final example, let us assume that the sLDF has a 
narrow peak in addition to the MS contribution. This scenario 
describes the case in which PMS stars have all the same lu- 
minosity, as in the horizontal branch or in the red giant clump 
phase; this case can be modeled by adding a pulse function to 
the MS section of the sLDF. Since PMS luminosities are larger 
than MS luminosities, the pulse function is located at { max . The 
sLDF is therefore: 



If 1 +y6 - a > 0, the mean luminosity is driven by ^ max . In a 
typical situation with £1 as 3, the most luminous stars will dom- 
inate the luminosity if a < 4: this is the case of Salpeter's IMF 



A 

l-a 



(Anax ^TO ^ 

-4 



if £ £ (Anin, (TO), 
J if ( G ((to, (mwd- 



(17) 



6 



M. Cervino, V. Luridiana, and N. Cervino-Luridiana: Probabilistic synthesis models 



Again, we are forcing the PMS to have a fixed weight with 
respect to the MS phase for simplicity reasons. In this case, 

assuming again that £ T £ » tj^ : 

A A — — 

A*l — ~, 7, ^TO "i (Anax — ^tq ) Anax- (18) 

1 + p - a I - a 

So, again, the mean depends on both £ro and £ max . 

These three examples illustrate the general fact that the 
mean luminosity depends strongly on the value of ( max . The 
dependence on € max (t) would be e ven stronger for higher -order 
moments of the distribution. As iGilfanov et all (2004) point 
out, the dependence of the results of synthesis models on the 
high-luminosity end of the sLDF is so strong that, without such 
a limit, synthesis models could not even be computed! As a fur- 
ther example of the relevance of the upper limit of the sLDF, we 
have shown in a previous paper that it also defines a lower limit 
for the luminosity of real clusters to be des cribed by standard 
synthesis models dCervino & L uridiana 2004). 

These examples also illustrate the following interesting 
facts: ( i) The mean value does not necessarily give information 
on the distribution of luminosities, to the extreme that there 
can be situations in which the probability to find a source in 
the region around the mean value is zero (Fig. [2 bottom panel) 
. This is the opposite of a Gaussian distribution, in which the 
mean is also the most probable value. ( ii) Different distributions 
can have the same mean value but different variances. In fact, 
this circumstance permits to use surface brightness fluctuations 
(SBF: the ratio between the variance and the mean value of the 
luminosity function) to break the age-metallicity degeneracy, 
which makes clusters with different ages and metallicities have 
the same p' v . ( Hi) The value of ju[ - cr is negative in our three ex- 
amples. This shows clearly that assuming, e.g., that p\ + lcr in- 
cludes ~ 68% of the elements of the distribution can be grossly 
mistaken in the case of non-gaussian distributions. As we will 
see later, this can also be the case with the distribution function 
of the integrated luminosity of an ensemble: this limits the use 
of goodness-of-fit methods based on the comparison with the 
mean, such as the test. 

Summing up, our schematic but realistic sLDFs show that 
knowledge of the mean and the variance does not provide a 
handle on t he problem if t he sha pe of the distribution is not 
known. As Gilfan ov et al.l (|2004) argue, the use of the mean 
value instead of the most probable one produces a systematic 
bias in the determination of integrated properties. 

4.2. Implications of mass binning 

Although Eq. ^| is expressed as an integral, synthesis codes 
approach it through numerical approximations. There are sev- 
eral reasons for this: first, the input data are available in tab- 
ular format rather than as analytical relations; e.g., given the 
enormous complexity of stellar evolution it is not possible to 
derive stellar properties as analytical functions of, say, initial 
mass and stellar age. Rather, stellar tracks are computed for a 
discrete and finite set of mass values, and their properties are 
conveniently interpolated a posteriori. Furthermore, calculated 
stellar tracks are generally expressed in tabular form, although 



there have been sporad ic attempts at expressing them in ana- 
lytical form fe.g. iTout et alJfl996T) . Second, even if analytical 
relations existed, their integration would plausibly require nu- 
merical methods. Third, the numerical precision of computers 
is limited. These circumstances force the actual calculations of 
synthesis models to be numerical rather than analytical: as a 
consequence, the mass variable must be binned. We will focus 
here on the implications of mass binning for Eg. 1101 
Introducing binning, Eq.^|can be expressed as: 

Jr»m Llp 
e(m; t) m (m) dm * V w t (19) 
m iow t-r^ 

where the index i identifies the mass bin and fj(t) is the (time- 
dependent) luminosity of the 2-th bin; the approximation holds 
only if the luminosity £,-(f) is indeed representative of the whole 
mass bin. The coefficient w, is computed by means of the ex- 
pression: 

Wi = I ipM(m)dm, (20) 

where m: ow and m" p are the lower and upper limits of the 2-th 
mass bin (the specific way in which these limits are defined 
varies from code to code). In the framework of deterministic 
synthesis models, w, is deterministically interpreted as the frac- 
tion of the total number of stars enclosed in the given bin. In 
the probabilistic approach, however, such number is not fixed. 
Each star is either born within a given mass bin i, with a prob- 
ability or outside it. When iV tot stars are selected, the num- 
ber of stars in each mass bin follows a binomial distribution. 
Furthermore, since all the bins share the same number N tot of 
stars, there is a finite covariance among bins: that is, the col- 
lection of the mutually covariating binomial distributions con- 
stitutes a multinomial distribution. According to this approach, 
Wi also represents the mean (as opposed to exact) contribution 
of each bin to the total number of stars. This interpretation is 
shared by statistical standard codes and Monte Carlo codes. 

4.2.1 . Distribution of stars in each bin 

The statistical strand of standard models was born with the goal 
of devising statistical tools to be applied to synthesis models. In 
this context, the binomial probability distribution of stars in in- 
dividual bins had been a pproximated by a Pois son distribution 
to simplify its handling (iCervino et alJl200 2b). As is known, 
the Poisson approximation is accurate only when the number 
of events N is large and the probability p is vanishingly small, 
in such a way that Np stays finite: we will show in the follow- 
ing that such approximation is not always accurate enough for 
our problem. 

A first point to consider is that speaking of a distribution of 
stars in bins only makes sense if a value of N to t is considered. 
In both cases considered, Poisson and binomial, the mean value 
of the number of stars in a bin is //,•(/!,-) = Mot X w,. Let us now 
illustrate the difference between the two distributions for the 
case of N\ox = 1. The probability for the source to fall in the 
i-th mass bin is, in the Poisson approximation, pf(l) = W; e~ w '. 
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In the binomial case such probability is p b (l) = w,-, which is 
the correct value. The difference is non-negligible for large w, 
values, which is typical for bins at small masses. 

Furthermore, when 7Y tot stars are considered, the Poisson 
approximation predicts the ratio Mj/MJ of the variance to the 
mean to be unity, whereas in Monte Carlo simulations this ratio 
sho ws a trend toward smal ler values for small bin masses (Fig. 
1 of lCervifio et"ai1 E)02b). revealing the underlying binomial 
distributions, in which: 



M 2 MotW,-(l - wd 



m: 



Mot Wi 



= 1 - Wj. 



(21) 



Note that, as expected, the discrepancy with respect to 1 
grows larger f or larger w, values - i .e., bins at lower masses. 
Unfortunatelv. lCervino et ai]j200 2b) failed to interpret this re- 
sult in term of a binomial distribution, and proposed to solve 
the discrepancy by reducing the size of the bins, while keeping 
the Poisson approximation. Of course, reducing the size of the 
bin (or, equivalently, wj) leads to a closer similarity between the 
two distributions, but this solution is not always viable due to 
the limited resolution of the input ingredients. A further point is 
that, because Mot is shared by all the bins, the numbers of stars 
in different bins covariate: while this effect is neglected when 
the bin distributions are modeled by independent Poisson dis- 
tributions, it arises naturally when a multinomial distribution is 
used. 

The difference between the Poissonian approximation and 
the multinomial description is particularly clear if we compute 
the variance of the pLDF Consider the expression for the vari- 
ance obtained assuming that the number of stars in each bin is 
described by a Poisson distribution, and that there is no corre- 
lation between bins (e.g. lCervino et all20 02b. among others): 



M 



(22) 



Now, consider the expression corresponding to a multinomial 
distribution, where, by definition, cov(n,-, nj) = -N tot w,- wf. 



M 2 = Yjfiiindij + 22^' C0V( "'''"; ) 

i i i*i 

= Motxf^w,^ -Yj^iQ 2 )- 

* i i ' 

-MotX [YjYjMjWiWj ■ 



(23) 



Comparing Eq. [22] to Eq. [23] it can be seen that the lat- 
ter contains two additional terms: the term N tot x Z;( w iA') 2 
arises as a result of dropping the Poisson approximation, and 
the term C t €j cov(«,, nj) represents the mutual depen- 

dence of bins, and hence arises as a direct consequence of bin- 
ning. 

These considerations show that the Poisson approximation, 
motivated by simplicity of handling, may break down in our 
problem (see also lLucvll2000l) . While the distinction between 



Ga 





Fig. 2. Top: isochrones in the fy vs. B-V color-magnitu de dia- 
gram for three different age values computed by Girardi et al. 
( 2002]) based on the evolutionary tracks by iMarigo & Girardi 
(2001). Bottom: same as above, for Wie- 



the Poisson approximation and the multinomial description is 
important in order to make sense of Monte Carlo simulations, 
it is fundamental in view of the discussion on our proposed 
probabilistic formalism (Sect.|6). 

4.3. Fast evolutionary phases 

In writing Eq.|9j we have tacitly assumed that the luminosity is 
a well-behaved function of the mass so that the isochrone is al- 
ways defined. However, this condition is in fact often violated, 
since a typical isochrone features shallow sections as well as 
peaks and discontinuities in the (I, m) plane. Shallow sections 
correspond to quiescent phases of stellar evolution, where evo- 
lution is slow (e.g. the MS); peaks correspond to faster phases 
(e.g. the asymptotic giant branch, AGB); and discontinuities 
correspond to abrupt transitions between phases (e.g. the onset 
of central helium burning in intermediate mass stars) or jumps 
in stellar behavior (the transition between WR and non-WR- 
type structures). In the following, peaks and discontinuities will 
generically be referred to as fast evolutionary phases. Some of 
these cases are illustrated by the comparison of Fig.Elwith Fig . 
[3] In Fig. [2] the isochrones computed bv iGirardi et alJ |2002) 
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based on the evolutionary tracks 

by|M 

arieo & Girardil d20Q lh 4 
are shown, at selected age values. Fig.[3]shows the relation be- 
tween the initial mass and the luminosity in the V and K bands 
for the same isochrones, focusing on fast evolutionary phases. 

If the mean value of the sLDF is determined through the nu- 
merical approximation of Eq.[23 fast evolutionary phases are 
difficult to handle, because a small difference in initial stellar 
mass yields a large difference in luminosity, so that the result of 
the numerical integration crucially depends on which luminos- 
ity is chosen to be representative of a given mass bin. In princi- 
ple, this difficulty can be dealt with by decreasing the width of 
the mass bin and choosing mass bins that do not go across dis- 
continuities; however, the available resolution in mass is gov- 
erned by the evolutionary tracks used by the code, and is gen- 
erally too low to resolve adequately such phases in synthesis 
models. 

This is a severe problem, since fast evolutionary phases are 
ubiquitous in post main-sequence evolution, and at certain fre- 
quencies they bear a major weight in the luminosity balance. 
Unfortunately, the way in which this problem is tackled is of- 
ten labeled a 'technical detail' of the computation and dis- 
missed as unimportant, and thus the papers describing evolu- 
tionary synthesis models do not generally make any reference 
to its solution - in spite of its difficulty and of the potentially 
disastrous consequences of incorrect assumptions. Here are a 
few examples of the ways in which this proble m has been ap- 
proac hed: a) in the Starburst99 synthesis code JLeitherer et alJ 
Il999l) . for certain metallicity values, an undocumented stellar 
track at 1.701 M n is a dde d to the tabulated track o f 1 .70 M Q by 
ISchaller et ail (1 19921) and lSchaerer etalJ dl993al) . to avoid the 
mass bins going across the discontinuity of the stellar models' 
behavior at such mass (C. Leitherer, D. Schaerer, & G. Meynet, 
private communication); b) to deal with the same problem, ad- 
ditional evolutionary tracks around the same mass range are 
used i n the computation of the i sochrones bv ICastellani et all 
J2003I) and ICariulo et all d2004 (S. De glTnnocenti, private 
communication) and Br ocato et al.l (119 99) (E. Brocato, private 
communication); c) to avoid the intrinsic discontinuity in the 
isochrones, the sam e mass is used twice in the isochrones by 
iGirardi et alJ J2002I) . namely at the end of the red giant branch 
and at the beginning of the horizontal branch (S. Bressan, pri- 
vate communication). 

As a final point, consider that fast evolutionary phases are 
those stages in which the evolution is rapid in any of the rele- 
vant luminosity bands, and not only the bolometric luminosity. 
Therefore, an isochrone that is adequately sampled in one band 
is not necessarily so in all the bands. 

4.4. Transient phases and fuzzy stellar behavior 

Traditional synthesis models rely on the possibility to express 
the luminosity as a function of mass (Eq. |5Jl, so that the inte- 
grals over luminosity are transformed into integrals over mass. 
However, there are stellar phases in which this is impossible to 
do, because the luminosity at a given age is not univocally de- 

4 These models are available at 

http : / /pleiadi . pd . astro . it/ 



termined by the mass. A few examples are: the thermal pulse 
phase in AGB stars, which is poorly resolved by stellar evolu- 
tionary models; supernova (SN) light curves, in which a star of 
fixed mass spans a large range of luminosity almost instanta- 
neously; variable stars; rotating stars, in which the luminosity 
emitted depends on both the rotation velocity and the inclina- 
tion angle. Partly because of this difficulty, these phenomena 
are generally not included in synthesis models; whereas ap- 
proaching the problem from the point of view of distribution 
functions permits to account for them by including them in the 
sLDF - provided they can be modeled quantitatively. This point 
will b e developed further in Sect. 17.21 See also iGilfanov et all 
( 2004) for an example on the inclusion of variability. 

Of course, in many cases of interest the modeling of the 
phenomenon may be a difficulty in itself. For example, mod- 
eling rotation is in itself problematic. Rotating stellar models 
are just beginning to appear on the market, and the distribution 
of velocities in a stellar population is largely unknown. Stellar 
rotation is not yet included in synthesis models primarily be- 
cause of these uncertainties, and not only because synthesis 
models are not flexible enough. But if stellar rotation or any 
other distributed stellar behavior is properly understood, or at 
least satisfactorily modeled, our method will permit to include 
it in synthesis modeling. 

4.5. Implementation of model atmospheres 

Usually, the available model atmospheres form a coarse grid 
in the (log g, log T e ff) plane, whereas isochrones are gener- 
ally continuous in the plane. In order to assign a model atmo- 
sphere to each isochrone location, one can either choose the 
nearest atmosphere of the grid, or interpolate between nearby 
atmospheres. Assigning the nearest atmosphere implies a fur- 
ther binning of data and may originate jumps in the results, 
although it is often assumed that these jumps cancel on aver- 
age when a whole population is considered. This problem will 
not be furt her discussed here; a more ex tensive discussion can 
be found in lCervino & Luridianald2005i) . 

In the next section, we will review the ways in which the 
basic task of synthesis codes is accomplished, given the diffi- 
culties outlined here. 

5. Computational algorithms 

This section outlines the basic ways in which the task of com- 
puting the luminosity is specifically performed by standard 
codes. The first two parts. l5.ll and l5.2l describe techniques im- 
plemented in standard models. The third part, 15.31 describes 
how the Monte Carlo method is put into practice. 

5.1. Isochrone synthesis 

The commonest technique to compute the luminosity of stellar 
populations is known as isochrone synthesis. Isochrone syn- 
thesis is the numerical integration of the product between the 
IMF and the isochrone approximated as in the rightmost part 
of Eq. ^] /Uj(f) = Yji w i^i( t )- The strongest assumption of 
this numerical approximation is that each fj(t) is representative 
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Fig. 3. Details of fast evolutionary phases in the V (solid line, left axis) and K (dotted-line, right axis) bands. Bottom panel: 
complete isochrones. Middle panels: blow-up of the mass axis, showing details of fast evolutionary phases at three different ages. 
Top panels: same as middle panels, with a more extreme blow up. The isochrone set is the same of Fig- EJ 



of all the stages included in it. Fast evolutionary phases evi- 
dently pose a problem in the integration, since a small interval 
in mass can correspond to very different luminosity values. In 
isochrone synthesis, the problem is solved increasing the num- 
ber of mass bins that map fast evolutionary phases. However, 
this strategy eventually requires assuming mass bins narrower 
than the precision of the published tracks and isochrones. It is 
interesting to note that o ne of the main adva ntages of isochrone 
synthesis mentioned bv lCharlot & BruzualHl99ll) is that it pro- 
duces 'smooth' results; however, this advantage only hides the 
problem of fast evolutionary phases, but does not solve it. 

5.2. The fuel consumption theorem 

As an alternative method, one can avoid using the expression 
of Eq.[23 an d perform instead the integration in the luminosity 
domain, that is integrate lif\i€): in this way, the fast evolu- 
tion of luminosity is automatically taken care of. This is the 
basis of the so-called fuel consumption theorem (FCT): to un- 
derstand fully the FC T and its assumptions, we refer the reader 
to the original paper (Renzini & Buzzoni 1986; Buzzoni 1989, 
see also Marast onll998l) . 



Note that isochrone synthesis and the FCT can be coupled: 
the computation of isochrones c an be done so as to fulfill the 
FCT requirements. We refer to Mar igo & Girardi| 11200 II) for 
an ex tended discussion on the subject (see also lBressan et alJ 



5.3. Monte Carlo methods 



The Monte Carlo method can be implemented in different 
ways, which take advantage of its potential to varying degrees, 
so they are not completely equivalent. To understand this point, 
consider how a sample of stars selected through a Monte Carlo 
mass assignment can be handled: first, the stars can be either 
followed individually throughout their evolution, or grouped in 
mass bins characterized by average properties; second, in ei- 
ther case an atmosphere, used to transform the bolometric lu- 
minosity into observable quantities, can be assigned by either 
choosing the nearest model in the available grid, or interpolat- 
ing between nearby atmospheres (see !4.5> . Therefore there is a 
total of four ways in which a Monte Carlo-selected stellar pop- 
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ulation 5 can be treated. In the literature, there are examples of 
at least three of these: 

(i) The stars obtained through the Monte Carlo selection are 
grouped in bins; to each bin, the nearest available atmo- 
sphere model is assigned, which depends on the age con- 
sidered. That is, the total luminosity is computed by means 
of the expression 



(24) 



where «, is the number of stars that have fallen into the 
i-th mass bin, and 1; is the luminosity of the atmosphere 
assigned to the bin; Eq. [33] is analogous to Eq. ^] with 
the only difference that the n, values vary from simulatio n 
to simulation. This strate gy is followed bv lBruzu al (2002); 
Bruz ual & CharloJ J2003I) (G. Bruzual, private communi- 
cation). This implementation of the MC method is not 
very time-consuming, so it is quite effective in obtaining 
an overview of the effect of sampling in different observ- 
ables, which is in fact the goal of the quoted papers. It 
has, however, two important disadvantages: the mass bin- 
ning hinders a proper sampling of the discontinuities in 
the isochrones, and the l{ assignation may either smooth 
out or spuriously amplify transient spectral features (see 
Sect.lO. 

(ii) The stars obtained through the Monte Carlo procedure are 
followed individually throughout their evolution, but the 
luminosities are assigned as in the previous method, that 
is choosing the nearest model in the grid. This approach, 
although similar to isochrone synthesis, is fundamentally 
different in that it avoids the additional rebinning implicit 
in the isochrone synthesis method. Mode ls of this kind 



( 1991); Cervino & Mas-Hesse ( 


1994); Kurthet al. (1999); 


ICervino et alJ 


(2000). These models work better than the 



previous at mapping discontinuities in the isochrones (i.e., 
the value of dt/dm is be tter evaluated: see, e.g. Fig. 4 of 
Mas- Hesse & Kunthll99ll) . but present the same disadvan- 
tages with respect to the assignment of (i. 
(iii) A further possibility is to perform Monte Carlo simulations 
over the luminosity function: the Monte Carlo-sampled 
stars are followed individually, and each is assigned a tai- 
lored atmosphere. Doing so requires performing interpola- 
tions in the £, grid, and permits reproducing the path in the 
HR diagram of individual stars. The individual luminos- 
ity values obtained in this way are eventually summed to 
yield the total luminosity. These models exploit the poten- 
tial of Monte Carlo method to its maximum and are the 
only ones able to map correctly the luminosity function 
without the necessity of binning. A sufficiently numerous 
set of Monte Carlo simulations of this kind provides di- 
rectly the distribution function of the ensemble luminos- 
ity. The bleeding edge of this method lies, of course, in 

s A further degree of freedom in Monte Carlo simulations is that 
they can be performed by either fixing the number of stars or the to- 
tal stellar mass. This aspect will not be discussed here, and we will 
consider only Monte Carlo simulations with a fixed number of stars. 



the interpolations techniques used both in the stellar evo- 
lution and in the atmosphere assignment. Exam ples of ap- 
plicati o ns of this method can be found i n ICantiello et alJ 
(120031) : ICervino & Valls-Gabaud J2003I) : iGirardil (12000). 
among others. 

6. Old tools for a new approach: the probabilistic 
formulation 

The conclusions arising from this brief overview of population 
synthesis can be summarized in the following points: 

- Deterministic standard models are based on a misunder- 
standing of the computed quantities. While it is generally 
claimed that they compute the integrated luminosity of a 
model cluster, they in fact compute the mean luminosity of 
the sLDF. 

- Statistical standard models (i.e., those that compute the 
mean and the variance of the luminosity distribution) give 
the correct interpretation, but they have a limited interpre- 
tative power. Knowledge of the mean and variance is not 
enough to characterize a distribution, much less to explain 
its shape in physical terms. 

- Suitably-done Monte Carlo simulations have the potential 
to bypass most of the problems of population synthesis aris- 
ing as a consequence of the statistical nature of the problem. 
However, they are extremely expensive in terms of CPU 
time and disk storage space, and require a considerable hu- 
man effort to be analyzed. Furthermore and most impor- 
tantly, the analysis of a set of Monte Carlo simulations per- 
formed without a grasp of the underlying statistics would 
be purely phenomenological, precluding the possibility to 
generalize the results to further simulations or real clusters. 

In this section, we will describe a probabilistic formalism 
that automatically frames the problem in its most natural in- 
terpretation, and opens the way to a deep understanding of the 
underlying physical problem. The first subsection is devoted to 
recall standard statistical concepts for those readers without a 
background in statistics and probability theory. Those who are 
already expert in this field can skip this part. 

A note on definitions is due here. Throughout this paper, 
we adopt the conventional (though often overlooked) distinc- 
tion between probability theory and statistics. A statistical for- 
mulation is one that seeks to intepret a sample of experimen- 
tal data in terms of an underlying distribution. A probabilis- 
tic formulation is one that assumes an underlying distribution 
and uses it to predict the resulting distribution of experimental 
data. The traditional Monte Carlo method has followed a sta- 
tistical approach: conclusions were drawn from the observation 
of a spread in the results of simulations. The present paper lays 
the foundations for a probabilistic approach, in that it seeks to 
give a formal description of the underlying distributions, and 
to make quantitative predictions based on them. Obviously, the 
two approaches complement each other and partially overlap; 
concepts like the one of distribution function belong to both 
probability theory and statistics, thus we will refer to them as 
probabilistic or statistical alike. When it comes to define our 
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method with respect to existing ones, however, we will sys- 
tematically draw a distinction and refer to it as a probabilistic 
formalism. 



6.1. Basic statistical concepts 

As we have seen in Sect.|2 the luminosity of stars can be char- 
acterized in terms of an underlying sLDF (Eq.[5}- The exploita- 
tion of this approach permits finding an effective solution to the 
problem of computing the integrated properties of stellar pop- 
ulations. The formalism presented here to this scope is not new 
in terms of theory of distributions and can be found in an y ad- 
vanced textbook of statistics (e.g. Kendall & Stuart 1977). For 
completeness, we briefly review here the relevant concepts. 

The properties of the sLDF as a distribution can be studied 
by computing its moments; the first moment is the mean, which 
we already encountered in Eq.|5] 



Jo 



€<pl(€)<U. 



(25) 



The general definition of the n-th moment of the sLDF is the 
following: 



Jr*eo 
o 



(26) 



If a = 0, we call it 'raw moment', while if a — p' x we call it 
'central moment'; in particular, the mean luminosity, which is 
the main output of synthesis models, is the raw moment of 1st 
order. 

In the following, we will adopt the notation p n to indicate 
central moments and the notation p' n to indicate raw moments. 
Let us write down explicitly the expressions for the second cen- 
tral moment (or variance) and the second raw moment of the 
sLDF, and their mutual relation: 



Jo 

Hi = {t-y.\fi PL {€)d{ = 
Jo 



.a. 



= n 2 ~P-i > 

and, analogously, for the third and fourth moments: 
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(30) 
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(32) 



Let us also introduce the characteristic function of the 
sLDF, <p(p), that is its Fourier transform: 



Ip) 
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(33) 



The characteristic function has the following properties: the 
coefficients of its Taylor expansion in ip are the raw moments 
of the distribution, while the coefficients of the Taylor expan- 
sion in ip of its logarithm In (/>l(p ) are the so-called cumulants 
of the distribution, k„: 



ln< 



Sp) = 2^'<r— r . 



(34) 



It follows that moments and cumulants are related by the ex- 
pression: 
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(35) 



In particular, the relations between the first four moments and 
cumulants are the following: 



K2 = P2 =H2~p.f, 

K 3 = p.3 = P 3 ~ ^P'iP 2 + 2/^f, 



K 4 



p A - 3p 2 =p' 4 - 4p[p' 3 - 3/4 + yip x p 2 - 6yu'i 



(36) 
(37) 
(38) 
(39) 



An important property of cumulants is that they are inde- 
pendent of the assumed origin of the distribution, except for 
K\\ they are also called sometimes "semi-invariants" due to this 
property. If the origin is taken at the mean of the distribution, 
Ki = 0. 

Finally, the skewness, j\, and the kurtosis, j2, of the dis- 
tribution are defined through ratios of the third and the fourth 
central moments respectively to appropriate powers of the vari- 
ance: 
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(40) 
(41) 



(Note that the definition of skewness and kurtosis may vary 
from author to author; alternative definitions can be found 
at http://mathworld.wolfram.com) These two quantities 
enclose information on the shape of the distribution: the skew- 
ness gives an idea of how asymmetric the distribution is, and it 
can be related to the difference between the value of the mean 
and the mode (the most probable value). The kurtosis is a mea- 
sure of peakedness, i.e. of the symmetric deformation of the 
distribution with respect to a Gaussian. In particular, Gaussian 
distributions have y\ — and y 2 = 0; flatter distributions have 
negative kurtosis values, while peaked distributions have posi- 
tive kurtosis values. 
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6.2. From stellar to population luminosity functions 



In the previous section we have characterized the properties of 
the luminosity function of individual stars. Let's see now how 
the luminosity function of an ensemble of M to t sources can be 
computed. 

6.2.1 . Obtaining the pLDF: exact solution 

As a general rule, the PDF of an ensemble of variables is ob- 
tained as the convolution of the PDFs of the individual vari- 
ables. For example, let <p x (x) be the PDF of a variable x and 
ifiyiy) the PDF of a variable y independent of x. The probability 
density of a variable u = x + y is given by the product of the 
probabilities of y? x (x) and <p y (y) summed over all the combina- 
tions of x and y such that u — x + y: 



f 

•J — a 



<p x (z) <fy(u -z)dz = <p x (x) ® ifiyiy), 



(42) 



which is the definition of convolution. In our case, we are as- 
suming that all the stars have luminosities distributed following 
the same distribution function, <fi]_{(), and that the stars are in- 
dependent on each other. Therefore, the PDF of an ensemble of 
N\ox stars is obtained by convolving <fi\A() with itself Mot times: 



<PlJ£) = <Pl({) ® ® ... ® MO ■ 



(43) 



Hence, if the sLDF is known, the pLDF of an ensemble 
of Mot stars can be computed by means of a convolution. The 
convolution is conceptually straightforward, but it poses severe 
numerical problems. The reason is the following. The convolu- 
tion must be performed linearly in luminosity and, in the gen- 
eral case, the dynamic range in luminosities spans eight orders 
of magnitude, from 1CT 2 L to 10 6 L . On the other hand, most 
programs that perform convolutions are based on Fourier trans- 
form routines that require a set of points ordered regularly on 
the x-axis; each time a convolution is performed the number of 
points on the x-axis is doubled. So, for a resolution of, say, 0.01 
L Q (necessary to resolve the low end of the luminosity function) 
10 8 points would be needed to define the luminosity function, 
and this number would be doubled each time a convolution is 
performed. Therefore, the points necessary to compute even a 
very undersampled population (a moderate number of convolu- 
tions) diverge rapidly, making it numerically unfeasible. We do 
not know any computational routine powerful enough to per- 
form this task, and any feedback from the community about 
this subject is highly welcome. 

An alternative solution is making the convolution loga- 
rithmically in the Fourier space (see next section). However, 
the numerical computation of the Fourier transform of a func- 
tion with an irregular sampling is also difficult and, again, we 
haven't find any routine to perform it satisfactorily. 

6.2.2. Scaling properties of the LDF 

Convolutions in the normal space are equivalent to products in 
the Fourier space, so that: 



<f>uJP) = 0lCp) Mm , (44) 
In L|ol (P) = Mot x In <p L (p) = Mot x J] K r ^L. (45) 



Hence, the cumulants of the luminosity distribution of an 
ensemble of Mot sources, K„, can be easily obtained from the 
cumulants of the sLDF, k„, through the simple scale relation: 



K„ = Mot x K n . 

Since k\ = [i' x , this also implies that: 

M\ = Mot x fj.[, 



(46) 



(47) 



where MJ is the mean value of the distribution that describes 
the luminosity of an ensemble of Mot sources: in other words, 
the mean luminosity obtained by synthesis models is scalable 
to a cluster of any size - including one with only 1 source! 

Thus, probabilistic reasoning confirms the intuitive expec- 
tation that the properties of an ensemble of Mot stars can be 
obtained by direct scaling of the properties of the sLDF. This 
result is fundamental for the interpretation of the output of syn- 
thesis models in terms of real clusters. However, it is clear from 
the above derivation that this simple scaling rule only applies to 
cumulants, not to moments; it is cumulants that hold the scale 
relations between the properties of the sLDF and the distribu- 
tion of the total luminosity of an ensemble. But, by virtue of 
Eas. 1361 - 131"* cumulants can be related to moments: k\ is equal 
to the first raw moment, and the following ks are equal to the 
central moments of the same order, which can in turn be ex- 
pressed in terms of the raw moments: therefore, to characterize 
a distribution we can refer either to the moments or to the cu- 
mulants. 

It is immediately seen from Eg s . 1401 and 14 1 1 that the shape 
of the distribution of the ensemble, when expressed in nor- 
mal form (through a transformation of the distribution func- 
tion to one with zero mean and unit variance: X. — ♦ x — 
(X, — M[)/ y/Mi), can be easily related to the shape of the 
sLDF: 



r. 



and 



1 



■ 7u 



F2 = 72, 

Mot 



(48) 



(49) 



where F; and F2 are the skewness and the kurtosis of the distri- 
bution of the ensemble. Note that, in agreement with the central 
limit theorem, Fj — > and F2 — > for large enough Mot values, 
i.e. the distribution tends to a Gaussian. 

Although the previous relations are useful to unveil the 
scale properties of LDFs, knowledge of the moments of a dis- 
tribution is useful but not sufficient to analyze it if its shape is 
unknown. For most application, one needs to know whether the 
distribution can be approximated by a Gaussian, and in case 
it is not, which its shape is. The following section will deal 
with the problem of characterizing a distribution by means of 
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its cumulants. The technique suggested can be used to solve 
two different kind of problems: on one hand, it can be used to 
generate a theoretical pLDF from a sLDF when the convolution 
is not feasible. On the other, it can be used to infer the pLDF of 
an observed population. 

6.2.3. Obtaining the pLDF: approximate solution 

Alternative solutions go through obtaining an approximate ex- 
pression for the pLDF. A quantitative characterization of the 
pLDF by means of its cumulants can be obtained by means 
of approximate expressions. To this aim, we suggest using the 
Edgeworth's series, which can be written schematically as: 



ip(x) = Z(x) 



■2> 



(50) 



where x is the normalized luminosity defined above, Z(x) is 
the Gaussian distribution function, and the terms f, are ob- 
tained by the Chebyshev- Hermite polynomials multipl ied by 
powers of the cumulants (Blinni kov & Moessnerlll998l) . This 
series is a true asymptotic series, i.e. the error is controlled 
whe n the series is truncated to a finite number of terms n. 
As lBlinnikov & Moessnerl (fl998) demonstrate, the error is on 
the same order of the last term of the sum, t„. When the er- 
ror is small, i.e. the approximation is satisfactory, the term 
= Z"=i t i measures the deviation of the LDF from gaus- 
sianity. 

These properties can be used to obtain an explicit descrip- 
tion of the distribution and to estimate its degree of gaussianity. 
The algorithm to be followed is described here and represented 
in Fig. E| 

(i) The range of interest in x (i.e. the normalized luminosity) 
must be defined. This is necessary because convergence is 
first reached at small x, and propagates outward as more 
terms are included in the truncated series. 

(ii) The maximum deviation 6 from gaussianity must be cho- 
sen, in order to discriminate between non-Gaussian and 
quasi-Gaussian behavior. 

(iii) The maximum discrepancy e admissible between the trun- 
cated series and the LDF must be chosen. 

(iv) A truncated expression is computed with the first n terms. 

(v) At this point, \t„\ provides an estimate of the error. If |?„|/| 1 + 
2„ I > e, the error is too large, i.e. the truncated series is not 
a good approximation to the LDF: a further term must be 
added and the process resumed at step (iv). This step might 
require computing further cumulants, as higher-order terms 
of the series include progressively higher-order cumulants. 

(vi) As the number of terms retained increases, \t n \ becomes 
smaller than e and the expression progressively approaches 
the LDF, until it eventually becomes an acceptable ap- 
proximation. At this point, if |S„| < 5 the pLDF is quasi- 
Gaussian; otherwise, it is strongly non-Gaussian, a fact that 
must be taken into account when the distribution is ana- 
lyzed. In either case, the approximated expression can be 
used. 



Choose relevant range 
[-x, x] 



Choose maximum deviation 
from gaussianity 8 



Choose maximum error 
in the approximation 6 



Compute Edgeworth' s 
expression with n terms 




n=n+l 



LDF is non-gaussian 



LDF is quasi-gaussian 



Edgeworth's expression 
approximates LDF 



Fig. 4. Algorithm to obtain an approximate expression for a 
pLDF based on the Edgeworth's series. 



The algorithm is summarized in the flux diagram of Fig. |4j 
which permits to find an explicit analytical expression for a 
pLDF of unknown shape but known cumulants. 

As an example, we give here the explicit expression of the 
Edgeworth's series truncated to include terms up to n — 2: 



^L M (X) 



V2^ 



( 1 + - T] (x 3 - 3x) + — T 2 (x 4 - 6x 2 + 3) • 
\ 6 24 



1 

72 



— r^(x 6 -15x 4 + 45x 2 -15) 



(51) 



In the top panel of Fig.[5]we represent the region where Eg. 1511 
satisfies the first test of Fig. @J that is the range of Fj and F 2 val- 
ues in which Eq . 15 1 1 approximates the pLDF with an accuracy 
of 10% or better in a given interval of normalized luminosity 
of x. The dependence on x is represented in Fig.|5]by different 
shades of gray. 
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Compute r. , r. 




Convolution or Monte Carlo 
simulation needed 



t(x)sZ(x)(l+S 2 ) 
(non-Gaussianity) 



<f(x)sZ(x) 
(quasi-Gaussianity) 



LDF can be expressed 
analytically in [-3 a 3o] 



Fig. 6. Characterization of a pLDF based on Edgeworth's ap- 
proximation to the second order and a Gaussianity tolerance 
interval of + 10%. 



Similarly, the bottom panel of Fig. [5] describes the region 
satisfying the second test, i.e. the region where the pLDF can 
be approximated by a Gaussian within a 10%. As expected, this 
region is centered around the point [Fi = 0,T2 = 0], where a 
true Gaussian would lie. Again, this depends on the range of 
x considered: the wider this range, the narrower the range of 
acceptable Ti and F 2 values. 

The extension of the dark-gray region can be used to define 
a simplified diagnostic test for luminosity functions (Fig. |6j. 
This test is based on the algorithm of Fig.@J but only the first 
four moments are used, and conventional figures of 10% are 
used to define whether Edgeworth's approximation is accept- 
able and the LDF is quasi-Gaussian. Since Y\ and F2 are relat- 
able to the skewness and kurtosis of the sLDF via N lot (Has. 1481 
and!49ll. this test can also be used to determine the minimum 
number of stars necessary to ensure Gaussianity in a given 
band. An example of this technique will be given in Sect. 17. II 

6.2.4. Computation of the variance of the pLDF 



We have seen in Sect. l4.2l that the Poisson distribution is only 
an approximation to the distribution of the number of stars in 
a bin, and that it is not a safe one. The correct alternative is 
to use the multinomial distribution, which has the additional 
advantage of keeping track of covariance effects across bins. 
However, the multinomial distribution by itself does not tell 
the whole story. Additionally, its handling is difficult. We will 
show here that the adoption of the probabilistic formulation 
permits bypassing this problem, because it provides the tools 
to compute the relevant quantities without the need of making 



assumptions on the distribution of stars in bins, except its ran- 
domness 6 . 

As an example, let's use our probabilistic formalism to 
compute the variance of the distribution. To do so, we need not 
make any assumption on the distribution of the number of stars 
in bins: we just need to take into account the scale relations of 
luminosity functions (Eq. I46> and the properties of cumulants 
(Eqs.ESHS): 

M 2 = K 2 = N tot X k 2 = N tot x /i 2 . (52) 
If we apply the approximation of Eq.^3 we obtain: 

m 2 = #«x|jT *V(o«tf-(jT e<pL($dt 

= Mot x|^Wi^-|^W;£ J J= 

= Mo, x ( 2 - Yj&iti? - Yj X lit > Wi w j\ (53) 

' i i i j*i ' 

which is the same as the expression of Eq.[23] This simple ex- 
ample shows the power and, at the same time, the simplicity 
of our probabilistic formalism. Note, however, that this result 
does not imply that the multinomial description is wrong: on 
the contrary, it is implicit in the probabilistic treatment. Our 
point here is that the probabilistic treatment is simpler and more 
powerful than an analysis explicitly based on the multinomial 
distribution. 

Finally, note that the value of the variance of the pLDF cus- 
tomarily assumed in the literature (Eq. I22> is biased, with the 
work by ICantiello et alJ J2003I) as the only exception we are 
aware of. Note that in fact, the variance obtained from Eq. |22| 
is the second raw moment of the pLDF, and not its second cen- 
tral moment. 

6.3. Technical problems in the computation of the 
pLDF 

Unfortunately, the determination of the LDF of an ensemble 
is hindered by a few technical problems. The first one con- 
cerns the stellar LDF: although some aspects of stellar popula- 
tions (like variability or transient events) would be more easily 
treated via the luminosity function as compared to the standard 
method, the modeling of the luminosity evolution during fast 
evolutionary phases poses the same problems as in standard 
codes. 

A second issue is the contribution of dead stars to the lu- 
minosity function. Those stars have no influence on the com- 
putation of the moments of the distribution, except for a "nor- 
malization factor" that depends on the age of the population. 
In terms of the shape of the luminosity distribution, dead stars 

6 This is by itself, of course, a strong physical assumption. If, for 
example, mass segregation affects star formation, then this assumption 
would be false. In that case, the whole problem should be restated, 
since the mechanism of star formation should be incorporated in the 
Star Formation History. 
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show up as a pulse function located at zero luminosity. This can 
be easily understood in the following terms: let us assume an 
ensemble with iV tot initial stars. At a given age there is a non- 
zero probability that all stars are dead; the pLDF, marginalized 
to such a condition, is a delta centered on zero. There is also a 
non-zero probability that all but 1 stars are dead: the pLDF cor- 
responding to this case is the sum of a pulse function and the 
stellar luminosity function, the pulse function having in this 
case a strength smaller than the delta. Similarly, all but 2 stars 
can be dead, and the corresponding pLDF is a pulse function 
plus the convolution of two stellar luminosity functions. In to- 
tal, there are iV tot + 1 marginalized cases, corresponding to N to t, 
Mot - 1, N tot -2, 2 and 1 dead stars. The LDF of the ensemble 
should include all these cases, each with a weight given by its 
relative probability, so that it will have a pulse contribution cen- 
tered on zero. Although this result seems to contradict the cen- 
tral limit theorem, which states that the asymptotic shape of the 
LDF of the ensemble is a Gaussian, this is not the case, since 
the relative pulse contribution becomes smaller and smaller as 
Mot increases. 



7. Applications of the probabilistic treatment 

7.1. Characterization of a sLDF by means of its 
cumulants 



In Sect. I4.ll we showed that knowledge of the mean and the 
variance is not enough to characterize a LDF. As a rule of 
thumb, at least the first four cumulants should be taken into 
consideration, which describe the mean, the dispersion, the 
asymmetry and the peakedness of the distribution. As a first ex- 
ample of application of the probabilistic treatment, we will de- 
rive a few properties of stellar population from the analysis of 
its cumulants. Fig.0shows the first four cumulants of the sLDF 
for different bands and ages, obtained from the isochrones by 
iMarigo & Girardil fcOOll) . We have assumed a ISafoeted Jl955t) 
IMF in the mass range 0.15 - 120 M Q normalized by the to- 
tal number of stars. For a normalization in mass rather than in 
number of stars, these results must be multiplied by the mean 
mass, (m) = 0.52 M . The figure illustrates several interesting 
facts: 

- The value of ^/yuf = Hzlltf is in general larger than 10, 
with the exception of blue bands. As has been pointed out 
in Sect. 16.2.41 most synthesis models compute the second 
raw moment yL instead of the variance p,i. The difference 
between the two moments is a term in pf, but since pi is in 
general much larger than pf, p' 2 is numerically close to pi. 
Although the variance computed by most synthesis codes 
is systematically biased, in numerical terms the value is 
nearly correct. However, it is safer to compute the variance 
properly. 

- As pointed out in previous works, the Gaussian regime is 
not as common in stellar populations as often assumed. 
This is particularly the case of infrared wavelengths, whose 
values of y\ and j2 imply, given the criterion of Fig. [6] 
that at least 10 7 stars are needed (roughly m 5 x 10 6 M in 
initial mass) in the analysis building block (i.e. individual 



pixels in SBF studies, or the whole slit in integrated data) 
to ensure quasi-gaussianity. Non-gaussianity is a problem 
when one wants to obtain confidence intervals in terms of 
cr: when the distribution is non-Gaussian, it cannot be as- 
sumed that ju'j + 1 cr contains m 68% of the distribution. This 
is not necessarily a problem for goodness-of-fit tests like 
X 1 , since the test also works for non-Gaussian distributions 
provided the deviation from gaussianity is not severe. For 
example, T\ should be used to have an approximate idea of 
which test would be the best to compare the observations 
with the models. 



7.2. Inclusion of uncertainties in the input parameters 

The input ingredients used in population synthesis are gener- 
ally assu med to be fully k nown but are in fact affected by uncer- 
tainties JCervino & Luridianall2005l) . These uncertainties may 
either reflect incomplete knowledge or an intrinsic spread in the 
features of the quantities considered. In either case, they can be 
included in the sLDF provided they can be modeled quantita- 
tively. 

As an example, let's suppose that the mass distribution 
function depends on a parameter 9 that is itself distributed, that 
is let's replace the univocally determined function <p m (m) with 
a parametric function if/(m; 9). The parameter 9 can be charac- 
terized by a probability distribution function such that: 



p(9)d9 = 1, 



(54) 



where the integration interval is the range where is defined. 
To apply Eq. |9] it is necessary to eliminate the parametric de- 
pendence. This is done by integrating iff(m; 9) over all possible 
9 values, weighting it by p(9): 



-I 



<p m (m)= ilf(m;9)p{9)d9, 



(55) 



The resulting distributio n, which is technically called a mixture 
of p(6) and if/(m; 9) dKendall & Stuardll977h . can be used in 
Eq.fJl 

Let's consider two hypothetical examples. In the first, as- 
sume that if/(8; m) oc mT 6 and that the power-law index 9 has a 
rectangular distribution centered on Salpeter's index a = 2.35: 







2.35 < +S, 
otherwise, 



(56) 



with S > 0. If we mix if/(9; m) with p{9) and integrate, we find: 



<p m (m) 



- fm 



m) p(9) d9 oc 



(m d - m- d ) _ 235 

oc m . 

25 In m 



(57) 



When 6 — > 0, the expression {m 6 - m~ s )/(2Skim) — > 1 and 
Salpeter's law is recovered. When 6 is non-negligible, the sLDF 
is distorted with respect to the simple Salpeter's case (Fig. [8). 
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Age [a] Age [a] 

Fig. 7 . Main parameter s of the luminosity fu nction in sever al photometric bands, obtained from the isochrones by 
Mari go & Girardil feoOl) and adopting the IMF by Salpeterl( ll955b in the mass range 0.15 - 120 M . 



As a further example, let's assume that tff(6; m) oc m as 
above, but that now p(9) is gaussian: 

J _ (B-2.35) 2 

m = — 7= e ^ • (58) 

cr g yln 

In this case, the 0-weighted IMF is given by tp m (m) oc 
m~ 235 exp(criln 2 m) /2 and, again, the IMF is distorted with re- 
spect to the limiting case given by Salpeter's law (Fig.|8|l. 

Following this example, any spread in the input ingredients 
can be included in the modeling and contribute to the shape 



of the pLDF In particular, the method outlined above can be 
used to include transient phases and fuzzy stellar behavior in 
the modeling (Sect. 14.4b . For this to be done, it is just neces- 
sary that the phenomenon considered can be described in terms 
of a parameter of known distribution function. Other uncertain- 
ties that can be incorporated in the modeling are those that re- 
flect our imperfect knowledge of the problem: for example, the 
sLDF could be mixed with a Gaussian distribution to mimic 
observational errors. 
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Fig. 8. Salpeter's IMF (solid line) compared to two examples of 
distributed IMFs: rectangular (dashed, 6 = 0.5) and Gaussian 
(dot-dashed, cr = 0.5). 



It is important to note that the effect on the pLDF of a 
spread in the input ingredients cannot be determined a priori, 
as it depends on how the distribution of the input ingredient af- 
fects the sLDF In particular, it cannot be established whether 
the inclusion of distributed ingredients will render convergence 
to gaussianity more or less rapid. It is to be expected that, the 
further is the sLDF from gaussianity, the slower will the con- 
vergence of the pLDF to gaussianity be. For example, if the 
distribution of the input ingredient is bimodal, bimodality will 
persist in the pLDF of small N to t, and a larger number of stars 
will be required for the pLDF to converge to a Gaussian. 

7.3. Comparison between Monte Carlo simulations, 
numerical convolutions and the Edgeworth's 
approximation 

We have shown in Sect. 16.2. ll that obtaining the pLDF through 
a convolution of the sLDF is not simple; however, in simple 
scenarios it is possible to solve it. In the following, we will con- 
sider a simple case to illustrate this point. Our aim in this ex- 
periment is twofold: first, we will derive an explicit expression 



for the pLDF at different iV tot 's to show how its shape is related 
to the generating sLDF and how it depends on iV tot . Second, 
we want to show that, through our method, we are able to re- 
produce the main features of Monte Carlo simulations for any 
value of Ntot- 

Let us assume a stellar luminosity function made up of 
three Gaussians, representing the dead stars, the MS, and the 
PMS respectively. The parameters of the three Gaussians are 
chosen so that the broad features of a set of Monte Carlo sim- 
ulations for one star are reproduced (upper panels of Fig. |5J 
the vertical scale is logarithmic in the left panel and linear in 
the right panel). In particular, the mean and the variance of the 
triple-Gaussian LDF are constrained to be the same as those 
of the Monte Carlo simulations. We also confirmed a posteri- 
ori that j\ and yi are very close, which is expected for similar 
distributions. 

With a numerical routine, we convolved this sLDF with it- 
self Ntot times. The resulting pLDFs for selected values of iV to t 
are shown in Fig. |9] (solid line). The features of these pLDFs 
can be understood qualitatively as follows: the characteristic 
function of this sLDF is: 

MP) = A^e-^lP 2 -^ + A Mse -7W 2 -'W+ 

+Ap MS e^ <r ™ s '' 2 ~ , ' W , (59) 

where Ad s , ^ms, and Apms are the weights of the Gaussians cor- 
responding to dead stars, the MS, and the PMS respectively; 

^ms, and ^pms their locations on the luminosity axis; and 
o"d s , cr M s, and ctpms the respective dispersions. The exponent 
of the Mot-th power of this characteristic function is a sum of 
real exponents in p 1 (i.e. Gaussian distributions) and imaginary 
exponents in p (i.e. translations of the corresponding distribu- 
tions). Hence the final function will be a sum of Gaussians lo- 
cated at different positions. 

The vertical sequence of panels shows how, increasing N m , 
the pLDF progressively becomes smoother and more sym- 
metric, approaching a Gaussian shape. In the same figure, 
Monte Carlo simulatio ns with corresp onding iV tot values are 
also shown JCervino & Valls-Gabaudll2003l) . and these coin- 
cide remarkably with the analytical pLDF. This is a conse- 
quence of having chosen a sLDF similar to the Monte Carlo 
distribution function for Not = 1 (which, in turn, maps the 
underlying sLDF), and shows the power of the method: large 
Monte Carlo simulations become redundant, in the sense of be- 
ing predictable, if one can characterize the sLDF and succeeds 
in convolving it. 

Fig.^|compares the pLDFs obtained through convolution 
of the sLDF to their Edgeworth's approximation to the second 
order (Eg. 15 11 1. Each panel corresponds to a different number 
of stars; the cluster in the upper panel, with iV tot = 1000, is the 
same one of the lower panel of Fig. [9] The Edgeworth's approx- 
imation improves visibly across the range of iV tot considered: 
this is expected, because the closer is the pLDF to a Gaussian, 
the better is it approximated by Edgeworth's expression. 

Finally, in Fig.^2 we show a blow-up of Fig.|5]that com- 
pares the Fs of the pLDF obtained through the three different 
methods (Monte Carlo simulations, numerical convolution, and 
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the Edgeworth's approximation). Several points deserve to be 
emphasized here: 

- As noted above, the Ti and F2 values of the Monte Carlo 
simulations nearly coincide with those of the convoluted 
pLDF. This is a remarkable proof of the power of describ- 
ing stellar populations in terms of luminosity distributions. 

- The dots corresponding to clusters of 10 3 stars fall within 
the shaded region in the top plot (approximation test) but 
not in the bottom one (Gaussianity test), implying that the 
pLDF can be approximated by the Edgeworth's function, 
but is still far from Gaussianity. This is also apparent from 
the bottom panel of Fig.|9j which shows that the shape of 
the pLDF (described both as a convolution and by Monte 
Carlo simulations) is markedly multimodal and asymmet- 
ric, hence far from Gaussian. 

- Finally, in this example Gaussianity is marginally reached 
only above iV tot = 6 • 10 3 stars, and even a cluster as large as 
A^ot = 2- 10 4 is Gaussian only within the {—2a, 2cr] interval. 

7.4. Criteria for assessing the significance of fits 

As has been pointed out several times, the main result of 
synthesis models is the mean value of the stellar luminos- 
ity function, which can be scaled to clusters of any size. We 
have also shown that the relative dispersion of the model re- 
sults (the ratio cr(X)/M[ - yfM^/M^) decreases when N tot in- 
creases. However, cr(X.) increases in absolute terms, since it 
is proportional to VMot, a fact that should be taken into ac- 
count in the comparison of theoretical data to observed clus- 
ters. Furthermore, since each monochromatic luminosity has its 
own cr(X), they should be weighted differently in fits. Finally, 
although some regions of the spectrum may have a quasi- 
Gaussian distribution, this will not happen in general with all 
the regions. Hence, not all the frequencies (either in synthetic 
or in observed spectra) are equivalent, or even suitable, to ob- 
tain the properties of the observed cluster. 

As an example, Fig. 1121 shows the first four cumulants of 
a region of the visual electromagnetic spectrum for the pLDF 
of a 1 Ga cluster with solar metallicit y, obta ined from the syn- 
thesis code sed@ 7 using a SalpeterJ dl955l) IMF in the mass 
range 0.1-120 M Q . The results are normaliz ed to mass. The 
isochrones used are from Gi rardi et alJ (120021) covering a mass 
range from 0. 15 to 100 M n and b ased on the solar models 
(Z=0.019) by iGirardi et alJ I d2000b and iBertelli et alJ d 1994 
that include overshooting and a simple synthetic treatmen t 
of the thermal pulses AGB phase dGirardi and Bertellilll998j) . 
The atmosphe re models are t a ken from the high resolution li- 
brary by Martins et al. (2005); Gonzalez Delga do etaD J2005) 



7 sed@ is a synthesis code included in the Spanish Virtual 
Observatory and the Violent Star Formation Legacy Tool project by 
means of the PGos3 tool. The code is written in ANSI C under GNU 
Public License and, currently, is managed by M. Cervino; use of the 
code and its results must be referred to solely by its documentation 
(Sed@ Reference Manual, in preparation), its WWW address and the 
citations in the headers of the output VO Tables. The code is currently 
accessible on-line at http://ov.inaoep.mx/ Its inclusion in the 
VO service grid is under way. 



based on PHOENIX faauschildt and Baronlfl999t lAllar et alJ 
2001 ). the ATLAS 9 odfnew lib rary dCastelli and Kuruc3 
2003) computed with SPE CTRUM dGrav and Corballvll994lh 
the ATLAS 9 library dKuruczl 1 199 lb c omputed with 
SYNTSPEC faubenv. Lanz. and .Teffervll 1995b . and TLUSTY 
dLanz and Huben\l2003b at TFe/Hl ="o"o dex. 

The figure shows the region around H6. As expected from 
previous plots at these ages (Fig.Q, the values of Ti and T2 are 
quite low. However, not all the wavelengths have the same sta- 
tistical significance. In particular, the H<S line shows the largest 
values of Ti and T2, so in undersampled clusters its profile will 
be difficult to fit. On the other hand, the profile of Ca n H + He 
is a quite robust result, with a low relative dispersion and Fi and 
T2 values close to the continuum level. Finally, the Ca 11 K line, 
with Y\ and F2 values similar to those of the continuum, has 
a high relative dispersion. Summarizing, fitting the theoretical 
models of Ca 11 lines, including their profiles, to observed data 
would yield more realistic results than fitting either the inten- 
sity, the equivalent width or the profile of H<S. Of course, these 
conclusions depend on the age and metallicity. 

8. Future applications of the probabilistic 
treatment 

In this section we will briefly discuss several potential exten- 
sions of the formalism that could have a strong impact in the 
analysis of stellar populations. Details on a a few examples will 
also be given. 

8.1. Forthcoming extensions of the formalism 

Since in this paper we have only considered the case of SSPs, 
maybe the most important pending issue is the extension to 
other scenarios of star formation. This topic will be covered 
in a forthcoming paper. 

A current limitation of the formalism is that it only deals 
with integrated properties that scale linearly with the number of 
stars in the cluster. In the future, the formalism will be extended 
to include the case of luminosity ratios. To solve this problem, 
in addition to the cumulants of the distribution function of the 
ensemble, it is also necessary to obtain the correlation function 
of the corresponding quantities. 

A further assumption of the formalism in its present stage is 
that the stellar population has a fixed number of stars Mot- It is 
our intention to extend the formalism by including the case of 
a collection of populations with varying number of stars. This 
extended formalism could be applied to several problems, such 
as the analysis of stellar populations in pixels, which requires 
computing the global distribution resulting from the distribu- 
tion of stars in each pixel and the distribution of numbers of 
stars across pixels; the distribution of luminosities in globular 
clusters; the estimation of the difference between luminosity 
profiles in galaxies inferred by means of a comparison with the 
mode and the mean respectively; the comparison of theoretical 
SBF with observed ones; and the comparison among different 
Monte Carlo simulations performed with a total fixed mass. A 
few of these prospective applications will be discussed further 
in the remainder of this section. 
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Finally, the formalism developed here for the treatment of 
luminosity functions can also be applied to stellar yields. All 
the statistical considerations that we have made about luminosi- 
ties can be easily translated into equivalent issues in the field of 
chemical evolution, although in that case the star formation his- 
tory would have a fundamental role, and the differential equa- 
tions that describe (deterministically) the chemical evolution 
would become stochastic differential equations, whose mean 
would coincide with the results obtained deterministically. As 
in the present case, comparing the mean against observations 
can produce a bias that depends strongly on the shape of the 
distribution. 

The following sections will give more details on a few ex- 
amples among these. 

8.2. Surface Brightness Fluctuations 

SBF observations from galaxies and globular clusters have 
been proposed as a test of evolutionary tracks and isochrones 
(e.g. ICantiello et al1l2003[ among others). This test is based 
on the comparison between the observed variance across pixels 
in the image of a galaxy and the variance expected on statisti- 
cal grounds. However, there are several inconsistencies in this 
method as is applied at present: 

- Observational SBF are the result of an average over an 
additional distribution, the one that defines the number of 
stars falling in a given pixel. The theoretical SBF formalism 
doesn't take this second distribution in consideration. 

- Each of the building-blocks of the distribution (the pixels) 
is representative of the integrated luminosity of an ensem- 
ble. Although the formalism, in terms of moments, can be 
scaled to pixels containing any number of stars, the number 
of stars in a given pixel and the number of pixels with the 
same number of stars should be known, in order to evalu- 
ate the error in the estimation of the SBF due to the finite 
sampling of the underlying distribution. 

- The theoretical SBF used to date are computed under the 
implicit assumption of a SSP, whereas the mode of star for- 
mation of a real galaxy can be much more complex. 

This inconsistencies could be overcome by extending the 
formalism so as to include the LDF of populations with varying 
number of sources. This subject will be discussed at length in 
a forthcoming paper. 

8.3. Putting constraints on the globular clusters' 
distribution 

In a similar way, the luminosity distributions of globular clus- 
ters in galactic halos could be compared to the corresponding 
distributions obtained theoretically, either in terms of moments 
or in terms of the explicit shape of the distribution. The possi- 
bility of obtaining higher-order moments, such as the skewness, 
helps constraining the distribution, which is necessary to test 
evolutionary tracks and isochrones. However, there are draw- 
backs similar to those of the previous case, with the only ex- 



ception of the assumption on the validity of a SSP, which in a 
galactic halo is probably verified. 

To pursue the goals sketched above, it is imperative to know 
the initial distribution of cluster masses. In the following, we 
will show qualitatively that our method can also contribute to- 
ward this goal, by disclosing a potential source of bias in the 
use of synthesis models for the determination of the mass dis- 
tribution of globular clusters. Note, however, that a quantitative 
conclusion cannot yet be reached. 

In Fig. we show the pLDFs in Xv for clu sters wit h the 

LMC metallicity obtained from the isochrones by Girardi et al. 

J — — ^3, 

(2002) using the Edgeworth's approximation. We have plot- 
ted the distrib utions that correspon d to the extremes of the age 
ranges usedbv lZhang & Falll ifl999). The pLDFs correspond to 
cluster masses of 10 4 M (solid line), 10 4 5 M Q (dashed line), 
and 10 5 M (dotted line). The mean values of the correspond- 
ing pLDFs are marked as vertical lines. 

The figure shows that, given the asymmetry of the stellar 
luminosity function, most observed clusters will have smaller 
luminosities than the mean, i.e. smaller luminosities than those 
predicted by a standard code. The straightforward implica- 
tion is that the mass of observed clusters inferred by means 
of a comparison with standard models is, in most cases, un- 
derestimated, with a bias larger for smaller clusters (see also 
iGilfanov et all2004l) and younger ages. 

This effect is highly relevant for young and undersampled 
clusters. In partic ular, assum ing that the age estimation ob- 
tained bv lZhang& FalU (Tl 999) is not biased, the figure clearly 
shows that there can be a systematic underestimation of the 
mass of younger clusters, an hence, an overestimation of the 
number of low-mass clusters. In the cluster mass range consid- 
ered, this effect is more relevant in the 10 4 to 10 4 5 M interval. 
This fact could change the shape of the di stribution of the initial 
cluster masses obtained by Zhang ^ Falll (fl999 ). particularly in 
the low mass range, making it shallower or even inverting the 
slope. However, this result is only a qualitative application of 
the new formalism, and the example above should not be taken 
literally. To establish a firm conclusion it is necessary to apply 
the probabilistic treatment also to the determination of ages. 

8.4. Tracing the sLDF with resolved stellar populations 

The individual stars in a resolved population can be used to 
trace the sLDF. This can be done by comparing the first four 
cumulants computed by means of current synthesis models to 
the corresponding observed quantities. Resolved populations 
have several advantages: 

- iVtot is a known quantity (i.e. N tot -l, except in the case of 
crowding problems, which would be treated separately). 

- The number of stars used to map the luminosity function 
and the possible luminosity bias are known, and can be in- 
corporated theoretically (e.g. by changing the value of t m i n 
in the computation of the moments). 

- The assumption of a SSP can be directly tested. Moments, 
and their associated sampling errors, can be estimated from 
simple star counts using the standard formulations for unbi- 
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ased estimations of the mean, variance, skewness and kur- 
tosis. 

- It is possible to perform a tailored analysis consistent with 
the used isochrones without outsiders; e.g. Blue Stragglers 
can be excluded from the estimation of the moments if they 
are not included in the theoretical luminosity function. 

- Finally, the covariance effects between different evolution- 
ary phases are included in a natural way in the probabilistic 
treatment, both theoretically and observationally. 



8.5. Testing the isochrones 

The method proposed in Sect. I8.4l can also be used to compare 
the predictions of different sets of isochrones, since it is also 
sensitive to the number of stars expected in each region of the 
isochrone (which is almost IMF-independent for PMS s tars). A 
similar method has been proposed by Wilso n & Hurlevl (l2003). 

To illustrate the point, we show in Fig.^]the first four cu- 
mulants of th e sLDF for different band s and ages using the 
isochrones by Marigo & Girardl fcOQlh. computed following 
FCT requirements explicitly, and by Girardi et al. (2000) and 
iBertelli et ail Jl994 . We have assumed a lSalpeterldl955l) IMF 
in the mass range 0.15 - 120 M Q normalized by the total num- 
ber of stars, and the results have been obtained by a direct inte- 
gration of the isochrones. Note that, although the mean values 
and the fiz/fif' ratios are similar across the isochrone sets, there 
are large differences in Ti and T2 at some ages. This tells us, 
without even knowing the luminosity distribution function, that 
there are strong differences between both sets of isochrones, 
which can be directly related with differences in the treatment 
of stellar evolution, e.g. differences in the lifetimes of different 
phases. 

8.6. Application to the Virtual Observatory 

The method described in this work can be implemented in the 
VO as an automatized tool for the analysis of observed data, 
since synthetic models can be given in terms of probability dis- 
tributions suitable to be used in data mining algorithms or in 
Bayesian analysis. To achieve this goal, an appropriate theo- 
retical data model is necessary; the definition of such model 
is a task that is currently being carried out by the Theory in- 
terest group of the International Virtual Observatory Alliance 
(http://www.ivoa.neti. In addition, the extension of this 
probabilistic formalism to distributions of luminosity ratios, 
which are used in diagnostic diagrams, would be an asset for 
the development of more robust VO analysis tools. 

9. Conclusions 

This paper considers a series of problems in population syn- 
thesis that arise as a consequence of the distributed nature of 
stellar populations, and develops a new probabilistic formalism 
that takes them into account. With this formalism it is possible 
to reproduce and explain the features of Monte Carlo simula- 
tions without the need of performing them. The new formalism 
has several advantages with respect to Monte Carlo simulations 



in terms of generality and reliability. Unlike Monte Carlo simu- 
lations, it is not affected by sampling errors in the estimation of 
the moments. Unfortunately, its exact application requires com- 
puting repeated convolutions, and we do not know any compu- 
tational tool that can perform this task efficiently. Although the 
formalism is complete for luminosities, it must still be extended 
to the case of ratios. A summary of our conclusions follows: 

- In the first part, we revise the current standard formalism 
and discuss the phenomena it fails to address and the coin- 
cidences and differences with our method. The main differ- 
ence between the standard approach and ours is that the for- 
mer interprets the results of synthesis models deterministi- 
cally, whereas our formalism accounts for their statistical 
distribution. We show how the resulting distribution func- 
tions can be characterized in terms of their moments and 
cumulants, and how the shape of such distributions can be 
obtained from data, information that is necessary to com- 
pute the confidence intervals of the models' predictions. 

- Standard synthesis models work by handling the sLDF and 
not the pLDF Nonetheless, we show that pLDFs are the 
convolution of sLDFs, so that synthesis models can be used 
to study stellar populations, either integrated or resolved, 
once the correct probabilistic formalism is included. 

- The cumulants of the distribution scale with the size of the 
stellar populations. We explain these scale relations and 
show how the cumulants can be obtained in terms of the 
distribution moments. The variance as computed by syn- 
thesis models (generally in the context of SBF) is biased, 
as it actually corresponds to the second raw moment of the 
distribution, which is not a scalable quantity. Fortunately, 
given the high asymmetry and the power-law nature of the 
luminosity distribution, the numerical difference between 
the second raw and central moments is almost negligible. It 
is advisable, however, to use the right formula. 

- The standard deviation <x cannot be used as the unit of 
confidence intervals unless the distribution is known to be 
quasi-Gaussian. A nearly Gaussian regime is reached only 
when the sample contains more than > 10 5 stars, with the 
precise limit depending on the spectral region among other 
factors. It is therefore mandatory to perform safety checks 
before using tests that assume normality, such as the^ 2 test. 
In particular, clusters with large skewness should be treated 
with extreme caution, since in those cases the bulk of the 
distributi on is several crs aw ay from the mean value. We 
refer to lGilfanov etaD J2004 for a more detailed analysis 
of this issue. 

- The customary assumption of a Poisson distribution in bins 
is, in fact, not accurate enough. A realistic solution, fully 
consistent with the underlying distributions, is the multi- 
nomial distribution. The multinomial distribution describes 
in a natural way the covariance effects introduced by the 
binning, i.e. the correlations between different bins. 

- We give a few guidelines for assessing the robustness of 
fits and show that not all the features of the electromag- 
netic spectrum are equally suitable for using in the fitting 
of theoretical to observational data. Although the relative 
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dispersion <t(£)/M[ tends to for increasing iV tot , cr(-C) in- 
creases. 

- For practical applications, we show how the results of syn- 
thesis models can be directly applied to the study of re- 
solved stellar populations in quite simple terms, i.e. obtain- 
ing the distribution moments from the observed stars. This 
formulation allows reliable comparisons of observed data 
with theoretical stellar models to be made. 

- Current synthesis models cannot be used for comparison 
with observed SBF, since an additional distribution func- 
tion must be included in the treatment, describing the dis- 
tribution of number of stars across pixels. Our formalism 
provides a natural way to do it. 

- We give a preliminary example of how the probabilistic for- 
malism can be applied to the distribution of globular clus- 
ters and to chemical evolution models. The latter would re- 
quire the inclusion of stochastic differential equations. 

- Finally, we derived implications of this study for the devel- 
opment of analysis tools in the VO. For example, we are 
implementing synthesis models that include the probabilis- 
tic formalism in a new VO tool called PGos3, under devel- 
opment at http : //ov . inaoep . mx The present formalism 
will be also used for the development of data mining and in 
Bayesian algorithms. 

As our understanding of stellar populations shifts, popula- 
tion synthesis tools evolve. The problem of predicting the in- 
tegrated properties of stellar populations was initially framed 
as a deterministic one and solved by standard codes. A grow- 
ing awareness of the spread in the input parameters has boosted 
the interest in Monte Carlo simulations, whose phenomenolog- 
ical exploration has brought about important insights into the 
statistical aspect of the problem. The time is ripe now for a 
further forward step, one that advances the problem from a sta- 
tistical to a probabilistic formulation. As this evolution takes 
place, however, it is important to keep in mind that the new 
formalism reinterprets previous conceptions, rather than over- 
throwing them, and that it does not supersede the old tools, but 
instead aims to specify how and when they can be applied. The 
probabilistic formalism is best seen as a unifying model that 
includes the old tools and empowers them, in a direction that is 
becoming imperative for understanding the new observational 
data. 
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Appendix A: Notation 

This section contains an exhaustive summary of the notation 
used and its rationale. 

A.1. Quantities related to stars 

U is the luminosity of an individual object in a discrete sum 
(Eq.Q). 

t is a continuous variable representing the possible luminos- 
ity values of individual objects. 
fL = <fh(fi',f) is the probability density function (PDF) of the 
variable (, i.e. the stellar luminosity distribution function 
(sLDF) (Eq. |2J|. It depends on the assumed age (f, which 
is explicited here as a parameter), but also on the metallic - 
ity, the evolutionary tracks adopted, and the star formation 
history. 

4>l = <Pl(p\ t) is the characteristic function of tpi, defined as its 
Fourier transform (Ea.l33>. 

a4 = i s tne n " tn raw moment of the sLDF (Eas. l27H32> . 

l*n = f*n(t) is the n-th central moment of the sLDF (Eqs. 
!28H32b . The second central moment is the variance of the 
sLDF, and its square root, <x, is the standard deviation of 
the probability distribution. Note that in the literature the 
symbol <x is often used to represent the standard deviation 
of the sample, or mean square error, which is a statistical 
quantity. In this paper, we do not use <x with this statistical 
meaning. 

K n = K„(t) is the n-th cumulant of the sLDF (Ea.l35>. 
7n = Jnif) are the skewness and the kurtosis (72) of the 
sLDF (Eqs.goHlTJ. 

A.2. Quantities related to the integrated properties of a 
stellar population 

these quantities depend on the total number of stars Ntot in 
population, as well as on the age of the population: 

is the integrated luminosity of a given individual population 
(Eq.EJ. 

= -£(?) is a continuous variable representing the possible 
values of the integrated luminosity of stellar populations 
with N tot stars (Eq.[5}. 

= (fiL tM (£.', t) is the probability density function of X, i.e. the 
population luminosity distribution function (pLDF). Under 
suitable assumptions, it can be obtained by convolving the 
sLDF M ot times with itself (Sect. l6.2.lT >. 
= (p^iP; t) is the characteristic function of <£>l,„,(-0> defined 
as its Fourier transform (Ea.l44>. 
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M' n - M' n (t) is the n-th raw moment of the pLDF. 

M n - M„(t) is the n-th central moment of the pLDF. If n — 2, 

its square root is the variance of the pLDF, denoted as cr(X.) 

in this paper. 

K„ = K n {t) is the n-th cumulant of the pLDF. Under suitable 
assumptions, simple scale relations hold between K„ and 
K n (Eq.|46j. 

r„ = r„(?) is the skewness (Fi) and the kurtosis (Tz) of the 
pLDF. Under suitable assumptions, simple scale relations 
hold between T„ and y„ (Eas. l48H49l 



A.3. Quantities related to the IMF 

<Pm = <PM(tn) is the IMF, defined as the probability density func- 
tion for a star of having an initial mass m (Eq.[8}. 

<p' M = <p' M (m) is a function proportional to (fM(m): <p' M (m) = 
<PM(m)/(m), where (m) is the mean mass of the IMF (Eq. 
O. 

Wf = w;(m] ow , m" p ; f) is the probability that a given star has 
a mass belonging to the interval [mJ. ow , m" p ] (Eal20t. Note 
that the interval is defined arbitrarily depending on the com- 
putational needs and the characteristic evolutionary time of 
the population considered, so in most cases it is varied de- 
pending on the age of the stellar population. 

n, = «,(mJ ow , m" p ; t) is a random (distributed) variable rep- 
resenting the number of stars in a given mass interval 
[m' ow , m" p ] (Eas. 12211241 . Its mean value is <n,-> = Mot X w t . 
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Fig. 9. Monte Carlo simulations of clusters at solar metal- 
Mcity and 5.5 Ma in the K ba nd for different values of Af to t 
dCervino & Valls-Gabaudl2003[ shaded histograms), compared 
to analytical pLDFs obtained convolving A^ tot times a sLDF 
made up of three Gaussians (solid line). The vertical scale is 
logarithmic in the left panels and linear in the right panels. The 
analytical sLDF has the same mean and variance of the Monte 
Carlo simulations; the position of the mean is shown by a ver- 
tical dashed vertical line. 
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Fig. 11. Fi and F2 values of the pLDFs obtained through 
Monte Carlos simulations (black diamonds), sLDF convo- 
lutions (empty circles), and the Edgeworth's approximation 
(empty squares). Larger T\ and F2 correspond to lower N tot val- 
ues. 



Fig. 10. Comparison between the pLDFs obtained through 
sLDF convolutions (dashed line) and the Edgeworth's approx- 
imation (solid line) for clusters with different N iol values. 
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Fig. 12. Cumulants of the monochromatic spectral en ergy dis- 
tributi on for a 1 Ga luminosity function assuming a Salpeter 
ill 9551) IMF in the mass range 0. 15 - 120 M . The figure shows 
the region around H<5, which is marked with a solid line. The 
lines Ca n H + He are marked with dashed lines, and the line 
Ca ii K with dotted lines. 
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Fig. 14. Main parameters of the luminosity function for sev- 
eral photometric bands obtained fro m the isochrones b y 
i Marigo & Girardil(l200ll) (bold lines) and lGirardi etaD EoOO): 
Bertelli e t alJ ( 1 1994 ("light lines), assuming a Salpetea ( fl955T) 
IMF in the mass range 0. 15 - 120 M . 
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Fig. 13. Luminosity distribution functions in V for clusters with 
the LMC metallicity, approximated by the Edgeworth's series, 
for several ages and cluster masses: 10 4 M Q (solid line), 10 4 5 
M (dashed line) and 10 5 M Q (dotted line). The mean values of 
the pLDFs are marked as vertical lines. 



